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ABSTRACT 

The dynamical behaviour of near by trajectories is being estimated by Lya- 
punov Characteristic Exponents(LCEs) in the Generahzed Photogravitational 
Chermnykh-Like problem. It is found that the trajectories of the Lagrangian 
point L4 move along the epicycloid path, and spirally depart from the vicinity 
of the point. The LCEs remain positive for all the cases and depend on the 
initial deviation vector as well as renormalization time step. It is noticed that 
the trajectories are chaotic in nature and the L4 is asymptotically stable. The 
effects of radiation pressure, oblateness and mass of the belt are also examined 
in the present model. 

Subject headings: Trajectory: Lagrangian Point: LCEs: Photograviational: Chermnykh- 
Like Problem:RTBP 



1. Introduction 



In present paper our aim is to obtain trajectories of L4 and is to estimate the rate of 
deviation for initially closely related trajectories in the modified restricted three body prob- 
lem model (as in Kushvah (2008, 2009a)) with radiation from Sun, oblateness of the second 
primary (massive body) and influence of the belt. It is supposed that the primary bodies and 
a belt are moving in a circular orbits about the common center of mass of both primaries. 
First time such problem was discussed by Chermnykh (1987) and its importance in astron- 
omy has been addressed by Jiang and Yeh (2004a). More generalized cases of the problem 
were studied by many scientists such as Jiang and Yeh (2004b), Papadakis (2004),Papadakis 
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(2005) and Jiang and Yeh (2006); Yeh and Jiang (2006). The effect of radiation pressure, 
Poynting-Robertson(P-R)drag and oblateness on the hnear stabihty and nonhnear stabihty 
of the L4(-5) have been discussed by Kushvah and Ishwar (2006) ; Kushvah et al. (2007a,b,c). 
In our article Kushvah (2010), we have described the design of the trajectory and analysis 
of the stability of coUinear point L2 in the Sun-Earth system. 

The first fundamental article about LCN's was written by Oseledec (1968) in their study 
of the ergodic theory of dynamical system and Benettin et al. (1980) presented explicit meth- 
ods for computing all LCEs of a dynamical system. Then Jefferys and Yi (1983) examined 
stability in the restricted problem of three bodies with Liapunov Characteristic number. 
First time Wolf et al. (1985) presented an algorithm with FORTRAN code that allows to 
estimate non-negative Lyapunov Exponents(LEs) from an experimental time series. San- 
dri (1996) have presented method for numerical calculation of Lyapunov Exponents for a 
smooth dynamical system with Mathematica[Wolfram (2003)] code. Tancredi et al. (2001) 
compared two different methods to compute Lyapunov Exponents (LEs). They have shown 
that since the errors are introduced in the renormalization procedure, it is natural to expect 
a dependency of the estimated LCEs with the number of renormalization performed in the 
sense that the smaller the step the worse the estimation. In his study they made conclusion 
that the two-particle method is not recommended to calculate LCEs in these cases where 
the solution can fall in a region of regular or quasi regular solution of the phase space. For 
a region of strong stochastically the LCEs calculated with the two-particle method gives 
acceptable value. 

This paper is organized as follows: In section 2, we state the model of the dynamical 
system and compute the trajectories of L4. Section 3 gives method to compute the LCEs, 
where subsection 3.1 presents the first order LCEs for various set values of parameters, time 
ranges and renormalization time steps. Section 4 presents comment about stability using 
trajectories of L4. Lastly, section 5 concludes the paper. 

2. Trajectory of L4 

It is supposed that the motion of an infinitesimal mass particle be influenced by the 
gravitational force from the two primaries(massive bodies) and a belt of mass Mi,. We also 
assume that infinitesimal mass does not influence the motion of the two massive bodies which 
move in circular orbit under their mutual gravitational attraction. Let us assume that m\ 
and 1712 be the masses of the bigger and smaller primary respectively, m be the mass of the 
infinitesimal body. The units are normalized by supposing that the sum of the masses to be 
unity, the distance between both massive bodies to be unity. The rotating frame normalized 
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to rotate with unit angular velocity and the time is normalized in such a way that the time for 
one period as a unit so that, the Gaussian constant of gravitational = 1. For the present 
model, perturbed mean motion n of the primaries is given by = 1 + + -^-j^^^, where 
T = a + b, a, b are flatness and core parameters respectively [as in Yeh and Jiang (2006)] 
which determine the density profile of the belt; where = (1 — + fi , A2 = -^7^ is the 

oblateness coefficient of 7712] r^, Vp are the equatorial and polar radii of 1712 respectively, r is 
the distance between primaries and the radius of the belt: a = — — is a mass parameter; 
gi = 1 — -pT- is a mass reduction factor and Fp is the solar radiation pressure force which is 
exactly apposite to the gravitational attraction force Fg. In a rotating reference frame the 
coordinates of mi and m2 are (— /x, 0) and (1 — /i, 0) respectively We consider the model 
proposed by Miyamoto and Nagai (1975), and equations of motion are given as in Kushvah 
(2008) and Kushvah (2009a): 



X — 2ny 



y + 2nx — D,y, 



(1) 
(2) 



where 



2 (1 - /x)gi(a; + /i) /i(a; + /i-l) 3 iJ,A2{x + /j, - 1) 
r? r3 2 To 



MbX 



_ Wi 

(r2 + 72)3/2 72- 



(x + u) ^, , . 

2 — + + yVf + X - ny 



1 



2 i^y ^i^My 

ny- 



2 4 



Mby 



2^2)3/2 J.2 



y 

—{{x + n)x + yy} + y + n{x + jj) 



n^{x'^ + y2) (1 - ij)qi 11 IJLA2 



+ 



+ — + 



+ 



I {x + u)x + yy 
+Wx { ^ ^ — — - n arctan 

2rf 



r2 ■ 2rl ' (r2 + T2)^/^ 

y 



X + fjL 



1 /^^(l ^l) 2 / I \2 I 2 2 / I i\2 I 2 

n = (x + //) + y , r2 = (x + // - 1) +y . 



Cd 



The parameter Wi is considered due to P-R drag[more review in Poynting (1903), Robert- 
son (1937), Chernikov (1970), Murray (1994) and Kushvah (2009b)]. Where ri, r2 are the 
distances of m from first and second primary respectively. The dimensionless velocity of the 
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light is supposed to be Cd — 299792458. Then from equations (1) and (2) energy integral is 
given as: 



where the quantity E is an energy integral related to the Jacobi's constant C(= —2E). 

For numerical computation of equilibrium points, we divide the orbital plane Oxy into 
three parts x < — /x, —/i < x < I — ^ and 1 — ii < x with respect to the primaries. 
For the simplicity, we set pL — 9.537 x 10~^,T = 0.01. The equilibrium points are given by 
substituting fl^ — Qy — 0^ and presented in figure 1 when gi = 0.75, A2 = 0.25, M;, = 0.25. In 
this figure the dark blue dotes present the position of ^4(5) : {x — 0.347988, y — ±0.70645), 
the light blue represent the coUinear equilibrium points Li : x — 0.753578, L2 : x — 1.14795 
and L3 : X = -0.788385 for which y = 0. 

The equations (1-2) with initial conditions x{0) = y(0) = ^, x'{0) = y'(0) = 

are used to determine the trajectories of L4 for different possible cases. At at time t — 0, 
the origin of coordinate axes is supposed at the equilibrium point. 

In the present model all the computed trajectories of the L4 follow approximately the 
same path described by an epitrochoid whose parametric equations are given as: 



where ai is radius of a fixed circle, bi is radius of rolling circle and di is distance form 
center of rolling circle to to the point (,t (t), y(t)) which forms a trajectory. It is evident from 
above equations that if di depends on time then orbit is unstable and trajectory moves 
spirally outward the vicinity of the initial point. 

When qi = 0.75, A2 = 0, the trajectory is shown in figure 2 with panels(a-d) for M5 = 0.0 
and panels(e-f) for = 0.0, where frames(a& d) < t < 50, (b&e) 50 < t < 100 and (c&f) 
< t < 200. It is clear from figure that if < t < 50 and M), — 0.0 the trajectory of L4 is 
similar to the curve described by epitrochoid (4, 5) for oi = 1/7, 61 = 1 = di. When t > 50 
then di becomes function of time t, and the trajectory moves spirally outward. When Mh — 
0.25, the trajectory follows the path correspond to parameters Oi — l,bi — -\/5(irrational) , 
di — 3. Here the value of hi is irrational number which shows that the motion is non periodic. 



E — - (^x'^ + y^) — Q,{x, y, x, y) — (Constant) 



(3) 






When gi = 0.75, = 0.25, figure 3 depicts the trajectory for L4 with frames (a-c) for 
A2 = 0.25 and frames(d-f) for A2 = 0.50. In frame(a) < t < 50, (h) 50 < t < 75 and 
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Fig. 1. — The position of equilibrium points. 




Fig. 2. — Trajectory of L4 when qi = 0.75,^42 = 0.0 in frame(a) blue curve for < t < 50, 
(b) red for 50 < t < 100 and (c) green for < t < 200. Frames(a-c) Mi, = 0.0, and frames 
(d-f) Mb = 0.25. 
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(c) < i < 77 while (d) < t < 8.2, (e) 8.2 < t < 8.3 and (f)0 < t < 9. It is clear from 
frames(a-c) that the trajectory moves along approximately epicycloid path, when t increases 
it departs form the vicinity of L4. The region of stability shrinks and trajectory moves along 
a single cusped epicycloid, then it departs far from the initial point. Hence oblateness effect 
is significant factor to reducing the stability region. 



3. Lyapunov Chciracteristic Exponents(LCEs) 



It is well known that, if LCE> for some initial conditions which indicates the trajectory 
of initial condition is unstable. If LCE= for some values of initial conditions the orbit is 
neutrally stable and which corresponds to regular motion. If LCE< 0, the corresponding 
orbit is asymptotically stable. Now suppose S* be a 4-dimcnsional phase space such that 
S = {X : X = [x{t) , y{t) , Px{t) , Py{t)]^^''^'^} , then the time evaluation of the orbit is governed 
by the equation 



X = f{X) 



dH dH 
dx dy 
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The dynamical system is described by the Hamiltonian H which depends on Jacobian 
constant and given by 



1 



H ^7,{pI+ pI) + n{yPx - xpy) - U (x, y) 



(8) 



where PxiPy are the momenta coordinates given by 



Px = 



dH 
dpx 



Py = 



dH 

dpy 



U{x,y) = n 



n^{xp- + y^) 



Consider v = {6x, 5y, 6px, Spy) be a deviation vector from initial condition X{0) such 
that 1 1 Vol I = 1. Then the variational equation is given 



v^Df{X).v 



(9) 
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where superscript t over partial derivatives of U indicates their respective values at t etc. 
Then the Lyapunov Characteristic Exponent is given by 



For numerical computation of LCEs we use method presented in Skokos and Gerlach 
(2010) and Skokos (2010). To avoid overflow in numerical computation, we partition the 
closed interval / = [to,Tmax] into rii sub intervals with time step At and the time to run 
from to Tmax i.e. 

P{I) = {0 = to, ti, t2, ^3, • • • tk-1, tk,---, tni = Tmax} , (12) 
then equation (11) can be written as 

ni 

X{v{t))= lim VlogQ;(ife), (13) 

nit— >oo ' 
fe=0 

where a{tk) — \\v{tk)\\. To determine first order LCEs in next section, we will use initial 

vector X(0) = (0.499046,0.866025,-0.866025,0.499046) for classical RTBP(gi = 1,^2 = 
0,Mb = 0) and X(0) = (0.337957,0.81415,-0.954676,0.39629) for modified RTBP( = 
0.75, A2 = 0.25, Mh = 0.25). As in figure 4, at each step v{tk) will be evaluated from (6,10) 
using X{tk-i) and unit deviation vector v{tk-i) = V(tfc-i)(say). 



3.1. First Order LCEs 

Now consider = LDi, = LD2, = LD^ and R^ = LD^ spaces such that 
LDi D LD2 D LDs D LD^. To find the first order LCE(Ai), (2 = 1,2,3,4), we choose 
initial unit deviation vectors from LDi\LD2: vn = (1/2, 1/2, 1/2, 1/2), V12 = (0, ^), 
"^13 = (0, 0, ^), f 14 = (0, 0, 0, 1). The values of LCEs arc presented in log-log plot figure 
5 for t = — 10000 when qi = 0.75, A2 = 0.25, = 0.25 left panel corresponding to 
At — 1 and right for At — 2. Initially the values of LCE(Ai) are different, they are shown by 
curves(I)-(IV) correspond to four vectors respectively, but when t increases they merge into 
a single curve. To obtain LCE(A2), we choose initial unit deviation vectors from LD2\LD^ 
such that V21 = (^, 4^, 0), V22 = (0, \, \, 0), V2^ = (0, 0, 1, 0). Figure 6 shows LCE(A2) 
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Fig. 3. — Trajectory of L4 when qi = 0.75, = 0.25, where frames(a-c) for A2 = 0.25 and 
frames(d-f) for A2 = 0.50. 





Fig. 4. — In Plot (I) only one step is used to obtain LCEs while in (II) more that one 
normalization steps are used and V(tk-i) denotes the unit deviation vector for all k. 



02468 02468 
logCT] log(T) 

Fig. 5.— LCE(Ai) when qi = 0.75,^2 = 0.25, = 0.25 and < t < 10000; curve (I) 

vn = (l/2,l/2,l/2,l/2),(II):t;i2 = (0,^,^,^), (III):t^i3 = (0,0,^,^) and (Wy.Vu = 
(0,0,0,1). 
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Fig. 6.— LCE(A2) when gi = 0.75,^2 = 0.25, = 0.25 and < t < 10000; curves (I) 
^21 = i^v TTV TR^ 0),(II):t;22 = (0, 0) and (III):t;23 = (0, 0, 1, 0). 




log(T) log(T) 



Fig. 7.— LCEs when < t < 10000, qi = 0.75,^2 = 0.25 and = 0.25; where LCE(A3):(I) 
V3i = (75, 75, 0, 0),(II):t;32 = (0, 1, 0, 0) and LCE(A4)(III):t;4i = (1, 0, 0, 0). 
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when qi = 0.75, A2 = 0.25, Mb = 0.25, where left panel corresponds to At = 1 and right for 
At = 2. 

Now for computation of LCE(A3), we choose initial unit deviation vectors from LD^\LD4^. 
The results are presented in figure 7 for qi = 0.75, A2 = 0.25, Mf, = 0.25 and < t < 10000 
with left frame for At = 1 and right for At = 2. In figures 8 and 9, curves are plotted 
when qi = 1,^2 = 0, = 0, where left panel corresponding to At = 0.1 and right for 
At = 1. In figure 8, curves are labeled as(I) V21 = (;^, 0),(II):i;22 = (O'TJ'Tl'O)' 

(III):f23 = (0,0,1,0) and in figure 9, curves are plotted for < t < 100, where (I) 
~ (72' "75' 0) '(-'--'■) -"^sa = (O5 1) 0) 0)- The curves are in wave form with decreasing ampli- 
tudes which tend to zero at infinity and curves become constant. 

To determine LCE(A4), we choose V41 = (1,0,0,0) from LD4 . The corresponding LCE 
is shown by curve (III) in figure 7:qi = 0.75, A2 = 0.25, Mf, = 0.25 with left frame for At = 1 
and right for At = 2. In figure 10, we consider qi = 1, A2 = 0.0, = 0.0 in which curve(I) 
represents renormalization time step At = 0.1 and (II) for At = 1. It can be seen that (I) 
is a smooth curve and (II) is initially stepped curve but both curves are initially increasing 
in nature and after certain time they become constant. The LCEs(Aj), (i = 1,2,3,4) are 
presented in Table 1 for initial point X(0) = (0.337957, 0.81415, -0.954676, 0.39629) and 
qi = 0.75, A2 = 0.25, Mb = 0.25. It is clear from figures and Table that all first order LCEs 
are positive for various set values of parameters and renormalization time steps. This shows 
that the present dynamical system is stochastic. It is also noticed that if Tmax is not very 
large the LCEs depend on the choice of the renormalization time step as well as initial 
deviation vectors while if Tmax is very large then LCEs depend on renormalization time 
step only. 
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Fig. 10.— LCE(A4) when qi = 1,^2 = 0.0, = 0.0 and v^i = (1,0,0,0) where curve(I) 
represents the renormahzation time step At = 0.1 and (II) for At = 1. 



Table 1. First order LCEs for initial point X{0) = (0.337957, 0.81415, -0.954676, 0.39629) when 

qi = 0.75, A2 = 0.25, = 0.25. 



logt logAl('Ull) logAl(lll2) logAl(lll3) logAl(lll4) logA2(l'2l) logA2(l'22) logA2(l'23) logA3(i;3l) logA3(ll32) logA4('y4l) 



At = 1 


log 10" 


0.403836 


0.446978 


-0.474575 


0.182934 


0.176488 


0.233361 


-1.33349 


0.420617 


0.561358 


-1.92815 


log 10^ 


1.78922 


1.79078 


1.31139 


1.75423 


1.72431 


1.71777 


1.70948 


1.8.3046 


1.86196 


1.53422 


log 10^ 


3.89211 


3.89209 


3.83842 


3.88798 


3.88433 


3.88314 


3.88489 


3.89796 


3.90179 


3.8656 


log 10^ 


6.17181 


6.17188 


6.16639 


6.17139 


6.17102 


6.17098 


6.1711 


6.17242 


6.17281 


6.16922 


log 10* 


8.47207 


8.47216 


8.47153 


8.47203 


8.47199 


8.47208 


8.472 


8.47213 


8.47217 


8.47188 


At = 2 


log 10° 


0.266298 


0.287879 


-0.594977 


0.195472 


0.101962 


0.124369 


0.00722875 


0.335259 


0.418618 


-0.849014 


log IQi 


1.17546 


1.17649 


0.719618 


1.1448 


1.11437 


1.10752 


1.11319 


1.21633 


1.24543 


0.936813 


log 192 


3.2087 


3.20867 


3.15533 


3.20463 


3.20098 


3.19979 


3.20168 


3.21453 


3.21832 


3.18239 


log 10^ 


5.47966 


5.47973 


5.47424 


5.47924 


5.47887 


5.47883 


5.47895 


5.48027 


5.48066 


5.47707 


log 10* 


7.77903 


7.77911 


7.77848 


7.77898 


7.77895 


7.77903 


7.77895 


7.77909 


7.77913 


7.77884 


At = 10 


log IQi 


0.030892 


0.030311 


-0.25883 


0.011177 


-0.0066539 


-0.013192 


-0.002721 


0.0587713 


0.076341 


-0.10142 


log 102 


1.67368 


1.67363 


1.6238 


1.66989 


1.66652 


1.66538 


1.66726 


1.67913 


1.68263 


1.64948 


log 10^ 


3.87816 


3.87824 


3.87278 


3.87775 


3.87738 


3.87734 


3.87746 


3.87877 


3.87915 


3.87559 


log 10" 


6.17039 


6.17047 


6.16985 


6.17035 


6.17031 


6.17039 


6.17032 


6.17045 


6.17049 


6.1702 


At = 100 


log 10'-^ 


-0.0374842 


-0.0374685 


-0.064786 


-0.039577 


-0.0414376 


-0.042026 


-0.0410299 


-0.0344616 


-0.0325176 


-0.0507718 


log 10^ 


1.66081 


1.66088 


1.65587 


1.66043 


1.66009 


1.66006 


1.66016 


1.66136 


1.66172 


1.65846 


log 10"* 


3.87675 


3.87684 


3.87621 


3.87671 


3.87667 


3.87676 


3.87668 


3.87681 


3.87685 


3.87657 


At = 1000 


log 10^ 


1.66081 


1.66088 


1.65587 


1.66043 


1.66009 


1.66006 


1.66016 


1.66136 


1.66172 


1.65846 


log 10* 


3.87675 


3.87684 


3.87621 


3.87671 


3.87667 


3.87676 


3.87668 


3.87681 


3.87685 


3.87657 


At = 10000 


log 10* 


-0.0452992 


-0.0452123 


-0.0455707 


-0.0453203 


-0.045339 


-0.0452492 


-0.0453349 


-0.0452687 


-0.045249 


-0.0453539 



Note. — 



The values of LCEs depend on time step for normalization more than the initial deviation vector. 
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4. Stability of L4 



Now we suppose that the coordinates {xi,yi) of L4 are initially perturbed by changing 



indicates the direction of the initial position vector in the local frame. For simplicity, it 
is supposed that e = 0.001 and = f . We solved (1, 2) numerically using above perturbed 
initial point and plotted figure 11 when A2 = 0.0, which shows that the orbit of test particle 
and its energy constant. When qi — 0.75 we have panels(I&II) and qi — 0.50 then (III&IV) 
in which (I&III) describe the trajectory and (II&IV) correspond to energy integral E. It 
is clear from the orbit that initially trajectory moves in epicycloid path described by (4, 5) 
without deviating far from L4 and energy constant remains negative; but after a certain time 
it moves spirally outward from the region and energy constant becomes positive. Here blue 
curves represent = 0.25 and red curves correspond to = 0.50. 

The effect of oblateness of the second primary is shown in figure 12 when Mb = 0.25, 
where (l&ll) correspond to qi = 0.75 and (III&IV) for qi — 0.50. Panels (I& III) show the 
trajectory of perturbed point L4 and (II&IV) describe the energy integral of that point. The 
blue curves correspond to A2 = 0.25 and red for A2 = 0.50. The trajectory of perturbed 
point follows the path described by epitrochoid (4,5), as time increases it moves spirally 
outward from the vicinity of L4. It is seen that the oblateness is a significant effect on 
the trajectory and the stability of L4. When A2 — 0.0 the L4 is asymptotically stable for 
the value of t which lies within a certain interval. But if oblate effect of second primary is 
present (A2 7^ 0), the stability region of L4 disappears for large values of A2. 



From relation A2 = we obtain y = ^5^2 + ^. This shows that if A2 increases 

means the ratio — increases consequently — increases. Then form Ryabov and Yankovsky 
(2006), it is found that as the attracting particle recedes that is the ratio ^ diminishes, the 
difference between the attraction of spheroid and that of a sphere will decrease and, if r(or ri) 
is very large in comparison with rg, the spheroid will exert a force that practically coincided 
with that of a sphere. If A2 = 0.5 i.e. very large value then rg > r(hypothetically) then 
second primary becomes a thin fiat disk. In this case the both primaries have no separate 
gravitational attraction so they act like a single body and its sphere of infiuence is common 
with very large radius, that attracts perturbed point. Since perturbed point L4 is supposed 
in the equatorial plane of second primary, so the attraction of the oblate spheroid (equatorial 
bulge) upon L4 at a given distance from the centre of primary is greater than that of a sphere 
of equal mass(A2 = 0.0) which has been proved by Moulton (1960). The effect of oblateness 
can be seen in frame(c) of figure 2 and in frame(I) of figure 12, where attraction of the 
equatorial bulge of the second primary increases with A2. Hence from above discussion we 
can say that if A2 — 0.50(hypothetically), the trajectory suddenly moves from the vicinity 



a;(0) = Xi + ecos(0),i/(O) = yi + esin(0) where = arctan 
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Fig. 11.— Stability of L4 when A2 = 0.0 with panels (I&II):gi = 0.75 and (III&IV):gi = 0.50 
in which blue curves correspond to Mh = 0.25 and red for = 0.50. 
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Fig. 12. — Effect of oblateness on the stability of L4 when = 0.25 panels (I&II)gi = 0.75 
(III&IV) qi = 0.50 in which blue solid curves for A2 = 0.25, red curves for A2 = 0.50. 
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of L4 as the time increase. 

5. Conclusion 

We have obtained the trajectories of L4 and its perturbed point, for various set values of 
parameters. It is found that the trajectories move along the epicycloid path upto a certain 
time then they move spirally outward from the vicinity of the point. Prom the first order 
Lyapunov Characteristic Exponents(LCEs), we have seen that the behaviours of trajectories 
are stochastic. It is also found that the radiation pressure, mass of the belt and oblateness 
are significant effects, they reduce the stability region and increase the stochasticity in the 
system. It is also found that if A2 = 0.50(hypothetically), the trajectory suddenly moves 
from the vicinity of L4 as the time increase. 
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